module SquareTriangularGridManipulation
    use TriangleElementManipulation
    implicit none


	integer, parameter :: numberOfBoundaries = 4

	integer, parameter :: prescribedTemperatureBoundary = 1
	integer, parameter :: prescribedHeatFluxBoundary = 2
	integer, parameter :: convectionBoundary = 3

	real(kind=8), parameter :: MaximumError = 1.0d-6

    integer, parameter :: MaxSteps = 100000

    type SquareTriangularGrid
        integer :: Nelx
        integer :: Nely

        real(kind=8) :: Lx
        real(kind=8) :: Ly

        real(kind=8) :: Basex
        real(kind=8) :: Basey

        type(TriangleElement), dimension(:), allocatable :: Elements
        real(kind=8), dimension(:,:), allocatable        :: Nodes
        integer, dimension(:,:), allocatable             :: Neighbours !< (Element, NeighbourNumber) obs: triangles have 3 neihgbours each

        integer, dimension(:,:), allocatable             :: BoundaryEdges
    end type

contains

    !> \brief Return the number of elements of the grid
    !!
    !! \param this - The triangular square grid of the
    !! \return n - the number of elements
    !!
	pure function numberOfElements( this ) result ( n )
		implicit none
		type(SquareTriangularGrid), intent(in) :: this
		integer :: n
		n = this%Nelx*this%Nely*2
	end function


    !> \brief Return the number of Boundary Elements of the grid
    !!
    !! \param this - The triangular square grid of the
    !! \return n - the number of elements on the boundary
    !!
	pure function numberOfBoundaryElements( this ) result (n)
		implicit none
		type(SquareTriangularGrid), intent(in) :: this
		n = 2*(this%Nelx+this%Nely)
	end function



    !> \brief Set the dimensional data for the grid.
    !!
    !! \param this - The Grid data structure
    !! \param Lx - The Length on the X Axis
    !! \param Ly - The Length on the Y Axis
    !! \param Nelx - The Number of elements over the X Axis
    !! \param Nely - The Number of elements over the Y Axis
    !! \param Basex - The Spatial Position of the Bottom-Leftest point on the X Axis
    !! \param Basey - The Spatial Position of the Bottom-Leftest point on the Y Axis
    !!
    subroutine SetGridData( this, Lx, Ly, Nelx, Nely, Basex, Basey)
        implicit none
        type(SquareGrid), intent(out) :: this
        real(kind=8), intent(in) :: Lx
        real(kind=8), intent(in) :: Ly
        integer, intent(in) :: Nelx
        integer, intent(in) :: Nely
        real(kind=8), intent(in) :: Basex
        real(kind=8), intent(in) :: Basey

        this%Nelx = Nelx
        this%Nely = Nely
        this%Lx = Lx
        this%Ly = Ly
        this%Basex = Basex
        this%Basey = Basey
    end subroutine



    subroutine GenerateStructuredGrid( this )
        implicit none
        type(SquareTriangularGrid), intent(inout) :: this

        call GenerateNodes(this)

        call GenerateElements(this)

        call GenerateNeighbours(this)
    end subroutine


    !> \brief Generate the Grid Points with a Rectangular-Regular aspect
    !!
    !! \param this - The Grid data Structure
    !!
    subroutine GenerateNodes( this )
        implicit none
        type(SquareTriangularGrid), intent(inout) :: this

        integer :: i, j
        integer :: statReturn

        real(kind=8) :: dx, dy

        allocate(this%Nodes( (this%Nelx+1)*(this%Nely+1), 2 ), stat=statReturn )
        if( statReturn /= 0 ) then
            print *, 'Unable to allocate memory for Nodes. NumbNodes = ', (this%Nelx+1)*(this%Nely+1)
            return
        end if

        dx = this%Lx/this%Nelx
        dy = this%Ly/this%Nely

        forall( i = 1:(this%Nely+1), j = 1:(this%Nelx+1) )
            this%Nodes( (i - 1)*(this%Nelx+1) + j ,1) = this%Basex + (j-1)*dx
            this%Nodes( (i - 1)*(this%Nelx+1) + j ,2) = this%Basey + (i-1)*dy
        end forall
    end subroutine

    !> \brief Generate the data structures for the elements of the region
    !!
    !! \param this - The Grid data structure
    !!
    subroutine GenerateElements(this)
        implicit none
        type(SquareGrid), intent(inout) :: this

        integer :: i, j
        integer :: statReturn
        integer :: elementCount
        integer :: square
		integer :: numberOfGhostElements

		numberOfGhostElements = this%Nelx + this%Nely + this%Nelx + this%Nely !< Number of elements around the perimeter

        allocate(this%Elements(numberOfElements(this)+numberOfGhostElements), stat=statReturn )
        if( statReturn /= 0 ) then
            print *, 'Unable to allocate memory for Elements. NumbElements = ', numberOfElements(this)
            return
        end if
		elementCount = 0
		i = 1
		j = 1
        do square = 1, this%Nely*this%Nelx !< Each square of the grid
			!< Setting triangle one
			elementCount = elementCount + 1
			call SetNodesTriangularElement(this%Elements(elementCount), &
						(i-1)*(this%Nelx) + j + 1, &  ! SE node
						(i)  *(this%Nelx) + j + 1, &  ! NW node
						(i-1)*(this%Nelx) + j )     ! SW node
            !this%Elements( elementCount )%Nodes(3) = (i)  *(this%Nelx) + j      ! NE node
            call CalculateNormalOuterVector
            !< Setting triangle two
			elementCount = elementCount + 1
			call SetNodesTriangularElement(this%Elements(elementCount), &
						 (i-1)*(this%Nelx) + j, &     ! SW node
						 (i)  *(this%Nelx) + j, &     ! NE node
						 (i)  *(this%Nelx) + j + 1 ) ! NW node
            !this%Elements( elementCount )%Nodes(2) = (i-1)*(this%Nelx) + j + 1  ! SE node
            call CalculateNormalOuterVector
            call CheckNextLine
        end do
	contains
		subroutine CheckNextLine !< Maximum number of elements in a line
			if ( j == this%Nelx ) then
				j = 1
				i = i+1
			end if
		end subroutine
		subroutine CalculateNormalOuterVector
			call CalculateOuterVector(this%Elements(elementCount), &
						this%Nodes(this%Elements(elementCount)%Nodes))
		end subroutine

    end subroutine


    !> \brief Generate the Neighbours for the Elements on the Grid. It does not include the Ghost Elements on the Boundary
    !!
    !! \param this - The Grid data structure
    !!
    subroutine GenerateNeighbours(this)
        implicit none
        type(SquareGrid), intent(inout) :: this

        integer :: i, j
        integer :: n
		integer :: square
        integer :: statReturn

        allocate(this%Neighbours(numberOfElements(this), nNodes), stat=statReturn )
        if( statReturn /= 0 ) then
            print *, 'Unable to allocate memory for Neighbours Elements. NumbNeighElements = ', numberOfElements(this)*nNodes
            return
        end if

        this%Neighbours = 0

        do i = 1, this%Nely
            do j = 1, this%Nelx
                n = (i-1)*(this%Nelx) + j
                if(i /= 1)         this%Neighbours(n,1) = (i-2)*(this%Nelx) + j      ! Southern Neighbour
                if(j /= this%Nelx) this%Neighbours(n,2) = (i-1)*(this%Nelx) + j + 1  ! Eastern Neighbour
                if(i /= this%Nely) this%Neighbours(n,3) = (i)*(this%Nelx) + j        ! Nothern Neighbour
                if(j /= 1)         this%Neighbours(n,4) = (i-1)*(this%Nelx) + j - 1  ! Western Neighbour
            end do
        end do
    end subroutine

!// \todo (ronaldo#1#): Gerar uma forma de alocar corretamente os Elementos Fantasmas
    subroutine GenerateBoundaryEdges(this)
        implicit none
        type(SquareGrid), intent(inout) :: this

        integer :: i, n

        integer :: statReturn

        n = 2*(this%Nelx+this%Nely)
        allocate(this%BoundaryEdges(n,2), stat=statReturn )
        if( statReturn /= 0 ) then
            print *, 'Unable to allocate memory for Boundary Edge Elements. NumbBoundEdges = ', (this%Nelx+this%Nely)*2
            return
        end if

        this%BoundaryEdges = 0

        do i = 1, this%Nelx
            this%BoundaryEdges(i,1) = i
            this%BoundaryEdges(i+this%Nely,1) = (this%Nelx*this%Nely) - (i-1)
        end do
        do i = 1, this%Nely
            this%BoundaryEdges(this%Nelx+i,1) = this%Nelx * i
            this%BoundaryEdges(i+2*this%Nelx+this%Nely ,1) = (this%Nelx*this%Nely) - i * this%Nelx + 1
        end do
    end subroutine

    !> \brief Return A Ghost Boundary Element Without any Degrees-of-Freedom setted
    !!
    !! \param boundaryElement - The Element which has the
    !! \param sharedEdge -
    !! \return Ghost - The Ghost Element Structure
	function GhostTriangle( boundaryElement, sharedEdge ) result (Ghost)
		implicit none
		type(TriangleElement), intent(in) :: boundaryElement
		integer, dimension(2), intent(in) :: sharedEdge

		integer, dimension(MatrixSize):: Indexes

		Ghost % Nodes(1) = 0
		Ghost % Nodes(2) = sharedEdge(1)
		Ghost % Nodes(3) = sharedEdge(2)

		Ghost % Matrix = 0.0d0
		Ghost % DoF = 0.0d0

		Ghost % NeighbourReNumbering = 0
		call SetNeighbourNumbering(Ghost, 1, boundaryElement)
		Ghost % S = 0.0d0
	end function


    !> \brief Set Data for Boundary 'Ghost' Triangle on the Heat Convection Boundary
    !!
    !! \param Ghost - The Boundary 'Ghost' Triangle for the heat convection
    !! \param BoundaryElement - The Triangle element on the boundary
    !! \param sharedEdge - The Edge which is shared with the elements
    !! \param convectionFactor - The convection heat coefficient
    !! \param EnvironmentTemperature - The temperature of the external environment
	subroutine SetGhostData_Convection( Ghost, BoundaryElement, sharedEdge, convectionFactor, EnvironmentTemperature )
		implicit none
		type(TriangleElement), intent(out) :: Ghost
		type(TriangleElement), intent(in) :: BoundaryElement
		integer, dimension(2) :: sharedEdge
		real(kind=8), intent(in) :: convectionFactor
		real(kind=8), intent(in) :: EnvironmentTemperature

		real(kind=8), dimension(2) :: T

		T(1) = BoundaryElement%DoF(Temperature_Initial+Ghost%NeighbourReNumbering(1,2)-1, Dof_Commit)
		T(2) = BoundaryElement%DoF(Temperature_Initial+Ghost%NeighbourReNumbering(1,3)-1, Dof_Commit)

		Ghost % DoF(HeatFluxX_Initial+1) = convectionFactor *(T(1) - EnvironmentTemperature) * Ghost%OuterVector(1,XAxis)
		Ghost % DoF(HeatFluxY_Initial+1) = convectionFactor *(T(1) - EnvironmentTemperature) * Ghost%OuterVector(1,YAxis)

		Ghost % DoF(HeatFluxX_Initial+2) = convectionFactor *(T(2) - EnvironmentTemperature) * Ghost%OuterVector(1,XAxis)
		Ghost % DoF(HeatFluxY_Initial+2) = convectionFactor *(T(2) - EnvironmentTemperature) * Ghost%OuterVector(1,YAxis)

		Ghost % DoF(Temperature_Initial+1) = T(1)
		Ghost % DoF(Temperature_Initial+2) = T(2)
	end subroutine


    !> \brief Set Data for Boundary 'Ghost' Triangle on the prescribed Dirichlet Boundary
    !!
    !! \param Ghost - The Boundary 'Ghost' Triangle
    !! \param BoundaryElement - The Triangle Element Adjacent to the Boundary
    !! \param sharedEdge - The Edge which they share
    !! \param Temperature - The Temperature on the Boundary.
	subroutine SetGhostData_PrescribedTemperature( Ghost, BoundaryElement, sharedEdge, Temperature )
		implicit none
		type(TriangleElement), intent(out) :: Ghost
		type(TriangleElement), intent(in) :: BoundaryElement
		integer, dimension(2) :: sharedEdge
		real(kind=8), dimension(2), intent(in) :: Temperature

		real(kind=8), dimension(2,probDim) :: Q

		Q(1, XAxis) = BoundaryElement%DoF(HeatFluxX_Initial+Ghost%NeighbourReNumbering(1,2)-1, Dof_Commit)
		Q(1, YAxis) = BoundaryElement%DoF(HeatFluxY_Initial+Ghost%NeighbourReNumbering(1,2)-1, Dof_Commit)

		Q(2, XAxis) = BoundaryElement%DoF(HeatFluxX_Initial+Ghost%NeighbourReNumbering(1,3)-1, Dof_Commit)
		Q(2, YAxis) = BoundaryElement%DoF(HeatFluxY_Initial+Ghost%NeighbourReNumbering(1,3)-1, Dof_Commit)

		Ghost % DoF(HeatFluxX_Initial+1) = Q(1,XAxis)
		Ghost % DoF(HeatFluxY_Initial+1) = Q(1,YAxis)

		Ghost % DoF(HeatFluxX_Initial+2) = Q(2,XAxis)
		Ghost % DoF(HeatFluxY_Initial+2) = Q(2,YAxis)

		Ghost % DoF(Temperature_Initial+1) = Temperature(1)
		Ghost % DoF(Temperature_Initial+2) = Temperature(2)
	end subroutine


    pure function MaximumError(this) result (rValue)
        implicit none
        type(SquareTriangularGrid), intent(in) :: this

        integer :: i, nEl
        real(kind=8) :: rValue

        nEl = numberOfElements(this)
        rValue = 0.0d0
        do i = 1, nEl
            rValue = max(rValue, SolutionErrorRange(this%Elements[i]))
        end do
    end function

	subroutine Step(this)
		implicit none
		type(SquareTriangularGrid), intent(out) :: this

		integer :: i
        integer :: stepsDone
        real(kind=8) :: maxerr

        maxerr = 1.0d6

        do while (maxerr > MaximumError .or. stepsDone <= MaxSteps)
            do i = 1, NumberOfElements(this)
                call StepElement(this%Elements(i), this%Elements(this%Neighbours(i,:)) )

            end do
        end do

	end subroutine

end module
